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We investigate aspects of universality of Glauber critical dynamics in two 
dimensions. We compute the critical exponent z and numerically corroborate 
its universality for three different models in the static Ising universality class 
and for five independent relaxation modes. We also present evidence for 
universality of amplitude ratios, which shows that, as far as dynamic behavior 
is concerned, each model in a given universality class is characterized by a 
single non-universal metric factor which determines the overall time scale. 
This paper also discusses in detail the variational and projection methods 
that are used to compute relaxation times with high accuracy. 
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I. INTRODUCTION 



Critical-point behavior is a manifestation of power-law divergences of the correlation 
length and the correlation time. The power laws that describe the divergence of the corre- 
lation length on approach of the critical point are expressed by means of critical exponents 
that are dependent on the direction of this approach, which may, e.g., be ordering- field-like 
or temperature-like. The exponents describing the singularities in thermodynamic quanti- 
ties can be expressed in terms of the same exponents. In addition to the exponents defining 
these power laws, another critical exponent, viz. the dynamic exponent z, is required for 
the singularities in the dynamics. This exponent z is defined by the relationship that holds 
between the correlation length ^ and the correlation time r, namely t cc 

One of the directions along which one can approach the critical singularity is the finite- 
size direction, i.e., one increases the system size L while keeping the independent thermo- 
dynamic variables at their infinite-system critical values. In this case, ^ oc L so that r oc L^. 
This relation has been used extensively to obtain the dynamic exponent z from finite-size 
calculations. 

In this paper we deal with universality of dynamic critical-point behavior. One would not 
expect systems in different static universality classes to have the same dynamic exponents, 
and even within the same static universality class, different dynamics may have different 
exponents. For instance, in the case of the Ising model, Kawasaki dynamics which satisfies 
a local conservation lawi has a larger value of z than Glauber dynamics,i in which such a 
conservation law is absent. Also the introduction of nonlocal spin updates, as realized e.g. 
in cluster algorithms, is known to lead to a different dynamic universal behavior. 

Conservation laws and nonlocal updates tend to have a large effect on the numerical 
value of the dynamic exponents, but until fairly recently, numerical resolution of the ex- 
pected differences of dynamic exponents of systems in different static universality classes 
for dynamics with local updates has been elusive. This is caused by the difficulty of ob- 
taining the required accuracy in estimates of the dynamic critical exponent. Under these 
circumstances it is evident that only a limited progress has been made with respect to the 
interesting questions regarding dynamic universality classes. 

In this paper we present a detailed exposition of a method of computing dynamic expo- 
nents with high accuracy.lifl We consider single spin-flip Glauber dynamics. This is defined 
by a Markov matrix, and computation of the correlation time is viewed here as an eigenvalue 
problem, since correlation times can be obtained from the subdominant eigenvalues of the 
Markov matrix. 

If a thermodynamic system is perturbed out of equilibrium, different thermodynamic 
quantities relax back at a different rates. More generally, there are infinitely many indepen- 
dent relaxation modes for a system in the thermodynamic limit. Let us label the models 
within a given universality class by means of k, and denote by the autocorrelation time 
of relaxation mode i of a system of linear dimension L. In this paper we present strong 
numerical evidence that, as indeed renormalization group theory suggests, at criticality the 
relaxation times have the following factorization property 

TLiK ~ m^AiL% (1) 

where is a non-universal metric factor, which differs for different representatives of the 
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same universality class as indicated; Ai is a universal amplitude which depends on the mode 
i; and z is the universal dynamical exponent introduced above. 

While the relaxation time of the slowest relaxation mode is obtained from the second- 
largest eigenvalue of the Markov matrix, lower-lying eigenvalues yield the relaxation times of 
faster modes. To compute these we construct, employing a Monte Carlo method, variational 
approximants for several eigenvectors. These approximants are called optimized trial vec- 
tors. The corresponding eigenvalues can then be estimated by evaluating with Monte Carlo 
techniques the overlap of these trial vectors and the corresponding matrix elements of the 
Markov matrix in the truncated basis spanned by these optimized trial vectors. It should 
be noted that both the optimization scheme and the evaluation of these matrix elements 
critically depend on the fact that the Markov matrix is sparse. That is, the number of 
configurations accessible from any given configuration is equal to the number of sites only, 
rather than the number of possible spin configurations. 

Given such fixed trial vectors, this approach has the advantage of simplicity and high 
statistical accuracy, but the disadvantage is that results are subject to systematic, variational 
errors, which only vanish in the ideal limit where the variational vectors become exact 
eigenvectors, or span an invariant subspace of the Markov matrix. Since the condition is 
rarely satisfied in cases of practical interest, a projection Monte Carlo method is then used, 
to reduce the systematic error, but this is at the expense of an increase of the statistical 
errors. The methods we use in this paper is a combination and generalization of the work 
of Umrigar et a/.,i and that of Ceperley and Bernu.i 

To summarize, the Monte Carlo method discussed here consists of two phases. In the 
first phase, trial vectors are optimized. The ultimate, yet unattainable goal of this phase is 
to construct exact eigenvectors. In this phase of the computation, very small Monte Carlo 
samples are used, consisting typically of no more than a few thousand spin configurations. In 
the second phase, one performs a standard Monte Carlo computation in which one reduces 
statistical errors by increasing the length of the computation rather than the quality of the 
variational approximation. 

The computed correlation times, derived from the partial solution of the eigenvalue 
problem as sketched above, are used in a finite-size analysis to compute the dynamic critical 
exponent z. We verify its universality for several models in the static universality class 
of the two-dimensional Ising model. We also address another manifestation of dynamic 
universality. As was mentioned, in addition to the usual static critical exponents, there is 
only one new exponent that governs the leading singularities of critical dynamics, viz.., z. 
Similarly, one would expect that, within the context of Glauber dynamics, the description of 
time-dependent critical-point amplitudes requires only a single non-universal metric factor 
to determine the time-scale of each different model within a given universality class. Our 
results corroborate this idea, which is the immediate generalization to critical dynamics of 
work on static critical phenomena by Privman and Fisher. 

In this paper we apply the techniques outlined above to three different two-dimensional 
Ising models subject to Glauber-like spin dynamics. These models are defined on a simple 
quadratic lattice of size L x L with periodic boundary conditions. The Hamiltonian Ti., 
defined on a general spin configuration S = (si, S2, . . .), is given by 



n{s) 



{ij) [ki] 



(2) 



kT 
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where the first summation is on all nearest-neighbor pairs of sites of the square Lx L lattice; 
the second summation is on all next-nearest-neighbor pairs; the Ising variables Sj, . . . , 
assume values ±1. Periodic boundaries are used throughout. In particular, we focus on 
models described by three ratios f3 = K'/K, namely f3 = —1/4, (nearest-neighbor model) 
and 1 (equivalent-neighbor model). The non-planar models for j3 ^ are not exactly 
solvable and their critical points are known only approximately. Yet, it was demonstrated 
to a high degree of numerical accuracy that that they belong to the static Ising universality 
class.til'EJ For the nearest neighbor model the critical coupling is i^' = | ln(l + a/2), for the 
other two models estimates of the critical points are K = 0.1901926807(2) for (3 = 1 and 
K = 0.6972207(2) for (3 = -l/4.Eili 

We use the dynamics of the heat-bath algorithm with random site selection. The single- 
spin-flip dynamics is determined by the Markov matrix P defined as follows. The element 
P{S', S) is the transition probability of going from configuration 5* to S'. If S and 5" differ 
by more than one spin, P{S', S) = 0. If both configurations differ by precisely one spin. 



P(y,5) = ^|l-tanh 



nis')-nis) 

2kT 



(3) 



where is the total number of spins. The diagonal elements P{S, S) follow from the 
conservation of probability 

En5',5) = l, (4) 

S' 

where S' runs over all possible 2^^ spin configurations. 

We denote the probability of finding spin configuration S at time t by pt{S). By design, 
the stationary state of the Markov process is the equilibrium distribution 

exp[-nis)/kT] _ Msr ... 

Poo[S) = = — - — , (5) 

where the normalization factor Z is the partition function. 

The dynamical process defined by Eq. (^ is constructed so as to satisfy detailed balance, 
which is equivalent to the statement that the matrix P with elements 

P{S',S)^^^^^P{S',S)MS) (6) 

is symmetric. Therefore the eigenvalues of P are real. 

The Markov matrix determines the time evolution of pt{S), i.e., 

Pt^,{S) = J2P{S,S')pt{S'). (7) 

S' 

The simultaneous probability distribution pt'^t'+t{S, S') that the system is in state S at time 
t' and in state 5" at time t' + t is 

p,,,+t{S,S')=P\S',S)ptiS), (8) 
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where P*{S', S) denotes the (5", 5") element of the t-th power of the matrix P. For sufficiently 
large times t', one may take pt'{S) = PooiS) so that the autocorrelation function Cyi(t) of 
an observable A, the average with respect to time t' of A{t')A{t + 1'), can equivalently be 
written as the ensemble average {A{t')A{t' + t)) for large t'. Thus 

CA{t) = lim Y: E A{S)A{S')pe,t'+t{S, S') (9) 

where ^(5*) denotes the value of A in a spin configuration S. After substitution of Eq. @ 
and expansion of A{S)pt'{S) in right-hand eigenvectors of the Markov matrix, it follows at 
once that the time-dependent correlation functions of a system of size L have the following 
form 

CAit) = ^csgnA^exph-^], (10) 

where the dependence on the specific model k has been suppressed in denoting by th 
the relaxation times of the independent modes of the equilibration process. The th are 
determined by the eigenvalues of the Markov matrix. We denote these eigenvalues Xu 
{i = 0,1,2, ... ,2^ — 1), and order them so that 1 = A^o > I-^liI > |-^L2| > ■ ■ ■• Note 
that conservation of probability implies that Xlo = by construction, the corresponding 
right-hand eigenvector is the Boltzmann distribution. 
The relaxation times are given by 

= = (11) 

The factor is inserted because, as usual, time is measured units of one fiip per spin, which 
corresponds to iterations of the process described by Eq. (0). 

Note that the stochastic matrix P has the same symmetry properties as the Hamiltonian 
and the Boltzmann distribution. In addition to spin inversion, these symmetries include 
translations, refiections and rotations of the L x L lattice. It follows that each eigenvector 
of P, as well as its associated relaxation mode, has distinct symmetry properties that can 
be characterized by a set of "quantum numbers." For instance, the eigenvector associated 
with the second-largest eigenvalue is antisymmetric under spin inversion, and invariant under 
translations, refiections and rotations. It describes the relaxation of the total magnetization; 
this process, with relaxation time r^i, is thus the slowest relaxation mode contained in the 
stochastic matrix. 

In this work, we restrict ourselves to relaxation modes that are invariant under geometric 
symmetries of the spin lattice. However, in addition to eigenvectors that are symmetric under 
spin inversion, we include antisymmetric ones, so as to obtain the longest relaxation time. 
As a consequence of this restriction to geometric invariance, spatially non-homogeneous 
relaxation processes fall outside the scope of this work. 

By design, the stationary state of the Markov process is the equilibrium state 

where the normalization factor Z is the partition function. 
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The dynamical process defined by Eq. (H) is constructed so as to satisfy detailed balance, 
which is equivalent to the statement that the matrix P with elements 



P(S',S)^^—PiS',S)MS) (13) 

is symmetric. 

The layout of the rest of this paper is as follows. In Section |I| we discuss the general prin- 
ciples of the method used in this paper. The expressions in this section feature exhaustive 
summation over all spin configurations, which renders them useless for practical computa- 
tions. The Monte Carlo summation methods employed instead are discussed in Section |T|. 



Trial vectors, a vital ingredient of the method, are discussed in Section |^ and numerical 
results are discussed in Section |V[ Finally, Section [V^ contains a discussion of issues that 
remain to be addressed by future work. 



II. EXACT SUMMATION 
A. Variational approximation 

1. A single eigenvector 

In this section we discuss trial vector optimization, which is used to reduce the statistical 
errors in the Monte Carlo computation of eigenvalues of the Markov matrix. First, we review 
the case of a single eigenvalue ,100 and then we generalize to optimization of multiple 
trial vectors. In this section, we discuss the exact expressions involving summation over 
all possible spin configurations. In cases of practical interest, these expressions cannot be 
evaluated as written; for their approximate evaluation one uses the Monte Carlo methods 



discussed in Section [II 



A powerful method of optimizing a single, many-parameter trial vector, say IV't), is min- 
imization of the variance of the configurational eigenvalue, which in the context of quantum 
Monte Carlo is called the local energy. That is, define iPt{S) = {S\i/jt) for an arbitrary 
configuration S. We wish to satisfy the eigenvalue equation 

^^(S) = A^t(^), (14) 

where the prime indicates matrix multiplication by P, i.e., f'{S) = Ss' -PI'S*, 5")/(S") for 
any function / defined on the spin configurations. Even if i/j^ is not an eigenvector, one can 
define the configurational eigenvalue by 

(15) 

I 0, otherwise. 

If ipT is not an eigenvector, Eq. (|^) gives an overdetermined set of equations for A, for 
a given ip^ and a sufficiently big set of configurations S. One can obtain a least-squares 
estimate of the eigenvalue A by minimizing the squared residual of Eq. ([T^) . This yields the 
usual variational estimate 
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which is the average of the configurational eigenvalue Ac- 

The standard Rayleigh-Ritz variational method, which can be used for the largest eigen- 
value consists in maximization of X{p) with respect to the parameters p. However, one can 
formulate a different optimization criterion as follows. The gradient of X{p) with respect to 
^t('S') is 

dX{p) _ ^ ^P'^{S) - X{p)MS) ^^^^ 



dMS) Es'MS') 

Clearly, this gradient vanishes for any eigenvector and this suggests as an alternative opti- 
mization criterion minimization of the magnitude of the gradient of a normalized trial vector 
ipT. With respect to Eq. (|T^), this corresponds to minimization of the normalized squared 
residual 

2( ^ Esms) - Mp)MS)? Esi^s) - x{p)?Msr {^t\{P - a)^|^t) . . 

which equals the variance of the configurational eigenvalue, as shown. 



B. Multiple eigenvectors 

Minimization of x^{p) is a valid criterion for any eigenvector, but if this is used with- 
out the equivalent of an orthogonalization procedure, one would in practice simply keep 
reproducing an approximation to the same eigenvector, the dominant one most of the time. 
Since orthogonalization is not easily implemented with Monte Carlo methods, we utilize 
straightforward generalizations of Eqs. (|1^) and (|16]) to deal with more than one eigenvalue 
and eigenvector. Eq. (|TS|) is a little problematic in this respect, as will become clear. 

Suppose we start from a set of n trial vectors ipTi, {i = 0, 1, ■ ■ ■ , n — 1). We can then 
write Eq. ( |14|) in matrix form 

n-l 

^^,(5) = ^A,,^T,(5). (19) 
i=o 

As before, the prime on the left-hand side indicates matrix multiplication by P, and again 
Eq. ([ToD for all i and S form an overdetermined set of equations for the matrix Aij. These 
equations have no solution, unless the n basis vectors ip^i span an invariant subspace of the 
matrix P, which in non-trivial applications of course is never the case. Again, however, one 
can solve for the matrix elements A^ in a least-squares sense. This yields 

A = VN-\ (20) 

where 

= Y.^T^iS)^PTAS) = (M^Tj) , (21) 

s 
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and 



h = E i^T^iS)^PTAS) = ii^Tmi^Tj). (22) 

5 

Note that although these matrix elements depend on the normalization of the \ipTi), the 
matrix A is invariant under an overall change of normalization. 

By diagonalization of the nxn matrix A one obtains an approximate, partial eigensystem 
of the Markov matrix. More specifically, suppose that 

A = D-Miag(Ao,...,A„_i)D. (23) 

The eigenvalues Aq > Ai > ■ ■ ■ > A„_i of A are variational lower bounds for the exact 
eigenvalues of the Markov matrix P, in the sense that A, < Aj,lll if the exact eigenvalues are 
numbered such that Aq > Ai > ■ ■ ■, in contrast with the convention used in the discussion 
following Eq. ([T0|). This property is a consequence of the interlacing property of the eigen- 
values of symmetric matrices and their submatrices, also known as the separation theorem.llZl 
Note that in denoting the eigenvalues we omit the index L indicating system size, where 
this is not confusing. The approximate eigenvectors tpi are given by 

n-l 

i'^ = Y. A,V'T„ (24) 

which can be verified as follows: multiply Eq. (|1^) through by D^i, and sum on i to verify that 
proportional to ipk- The expressions derived above are usuallyil derived by starting from 
the linear combinations given in this last equation. The Dij then are treated as variational 
parameters, and are determined by requiring stationarity of the Rayleigh quotient. This 
yields the following equation for the Dij 

^A,P,fc = A.^AAfc, (25) 
j j 

a generalized eigenvalue problem equivalent to the eigenvalue problem defined by A defined 
in Eq. 

Next we discuss the generalization to more than one trial vector of minimization of the 
variance as given by Eq. ([T8|) . In this context it is important to keep in mind that the vari- 
ational approximation is invariant under replacement of the basis vectors by a non-singular 
linear superposition. This yields a similarity transformation of A, and leaves invariant the 
approximate eigenvalues and eigenvectors. Of course, one would like to have an optimization 
criterion that shares this invariance. The squared residual of Eq. ([I9| ) fails in this respect. 
This sum is not even invariant under a simple rescaling of the basis functions t/'ti, and there 
is no obvious normalization comparable to the one used in Eq. (pISl). 

One way to perform the optimization in an invariant way is for each choice of the op- 
timization parameters to compute linear combinations I^j^jlV'Tj); where V is the nxn 
matrix, such that VAV'-^ is diagonal. Each of these linear combinations defines a via 
Eq. (|18|), and the parameters in the basis functions can then be optimized by minimiza- 
tion of these n sums of squares. One may define a convex sum of these xl and optimize 
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all parameters for all basis functions simultaneously with respect to this combined object 
function, or, as we did in our computations, one can perform the optimization iteratively 
one vector at a time for eigenvalues with increasing distance from the top of the spectrum. 

Another approach that also yields invariant results is to perform the optimization by 
dividing the set of configurations into several subsets, and computing a matrix A for each 
subset. One can then minimize the variance of the eigenvalues over these subsets. This is 
the procedure we followed to produce the results reported in this paper. 

We have not investigated which of the two procedures described above is superior. Both 
do have a problem in common, namely that for a wide class of variational basis vectors, 
they give rise to a singular or nearly singular optimization problem. This is a consequence 
of the fact that the basis states are not unique, even if the eigenvalue problem has a unique 
solution. For the optimization problem this means that there are many almost equivalent 
solutions, a problem commonly encountered in (non-linear) when on performs least-squares 
parameter fits. 

More specifically, if the basis vectors are such that a linear combination of trial vectors 
can be expressed exactly (or to good approximation) in the same functional form as the trial 
vectors themselves, then there is a gauge symmetry (or an approximate gauge symmetry) 
that yields a class of equivalent (or almost equivalent) solutions of the minimization problem. 
That is, if \ipTi) is a solution, J2j ^jlV'Tj) is an equivalent (or almost equivalent) solution for 
any V. This problem can be solved straightforwardly by fixing the gauge and performing 
the optimization subject to the constraint that J2j Vijli^Tj) oc \i'Ti)- If the gauge symmetry 
holds only approximately, this additional constraint may produce a sub-optimal solution. 



C. Beyond the variational approximation 

The eigenvalues obtained by the variational scheme discussed in the previous sections 
have a bias caused by admixture of eigenvectors in that part of the spectrum that is being 
ignored. This variational bias can be reduced in principle arbitrarily as follows.! 

Let us introduce generalized matrices with elements 

N,,{t) = {^T^\P'\^T,) (26) 

and 

A,(t) = (V^T.|^*+'|^T,). (27) 

For t = these expressions reduce to Eqs. (pi]) and (p2D. One can view the matrix ele- 
ments for t > as having been obtained by the substitution lipTi) P^^'^\4'Ti)- Expansion 
in the exact eigenvectors immediately shows that the spectral weights are reduced of "un- 
desirable" eigenvectors with less dominant eigenvalues, so that the vectors P^^'^\tpTi) span a 
more nearly invariant subspace of P than the original states. This process, however becomes 
numerically unstable as t — >■ oo, since in that case all basis vectors of the same symmetry 
collapse onto the corresponding dominant state. 
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III. MONTE CARLO SUMMATION 



Obviously, the summation over all spin configurations used in the expressions in the 
previous section can, in general, be done only for small systems. In this section, we discuss 
the Monte Carlo estimators of the expressions presented above. In principle, matrix multi- 
plication involves summation over all configurations and therefore is not practically feasible. 
However, for the dynamics we consider in this paper the summation required for the matrix 
multiplication by P in {S\P\'ipi) is an exception, since for a given S there are only con- 
figurations S' from which S can be reached with one or fewer spin flips, and these are the 
only configurations for which P{S, S') does not vanish. For all other configuration sums a 
Monte Carlo method is used. 

To produce a Monte Carlo estimate of ^is given in Eq. (|l^), sample spin con- 

figurations Sa with a = 1, . . . , Mc from the Boltzmann distribution ipBiSaY- This yields a 
Monte Carlo estimate of X{p) 



where 



MS.) = (29) 



Similarly, 



2. . ^ Eams.)-xMSa)? .3^. 



Parameter optimization for a single vector is done by generating a sample of a few thou- 
sand configurations and subsequently varying the parameters p while keeping this sample 
fixed. The same applies to the optimization of more than one vector, in which case estimates 
of the required matrix elements Nij and Vij are computed by 

-Y.^T^{SM^J{S^) = Ni„ (32) 

a 

and 

'P^J ~ E ^Ti{Sa)^Tj{Sa) = Vi,. (33) 
a 

We attached tildes to the symbols on the right-hand side of Eqs. ( |32| ) and (pSf ) to indicate 
that the corresponding quantities are stochastic variables, which is important to keep in mind 
for the following discussion. 

Since the matrix V is symmetric, one might be inclined to symmetrize its estimator Vij 
with respect to i and j. This symmetrization, however, destroys the zero- variance principle 
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satisfied by the expressions as written. As mentioned before, tlie eigensystem of A is obtained 
exactly and without statistical noise, if the basis vectors ip^i are linear combinations of n 
exact eigenvectors. In that ideal case, A is uniquely determined by Eq. ([I9|) even if it is 
applied only to an subset of configurations S. The same holds for a weighted subset as 
represented by a Monte Carlo sample. Even though the matrices V and N themselves 
depend on the weights and the subset, factors responsible for statistical noise cancel in the 
product VN~^. To demonstrate this, we write the estimator in matrix form 

N = M\ (34) 

and 

V = (35) 

where is a rectangular matrix with elements = 'ipTi{Sa), and ^-q, = 'ip'r^^{Sa)- Eq. (p!9D 
in matrix form becomes 

= A§. (36) 

Clearly, if this last equation holds, VN~^ = ^'^t(^^^t^-i = yY without statistical noise, as 
announced. 

We have assumed that one matrix multiplication by the Markov matrix can be done 
exactly; repeated multiplications rapidly become intractable. This is a problem for the 
computation of the matrix elements given in Eqs. ( p6D and (0). To obtain a statistical 
estimate of these matrix elements, one generates a time series with the Markov matrix P. 
One then exploits the fact that in the steady state of the Markov process, the relative 
probability of finding configurations Si, S2, ■ ■ ■ , St+i in immediate succession is given by 

P{St+i\St)---P{S2\Si)MSiY. (37) 
For a Monte Carlo run of length M^, this property allows us to write 

N^= E i^TiiSt)PiSt\St_i)---P{Si\So)i^TASo) 

So,---,St 

= E ^Ti{St)^TASo)P{St\St-i)---P{Si\So)MSof (38) 

So,---,St 

^ M-i E ^T^{S„+t)^Tj{S.) (39) 



cr=l 



Similarly, 



= E ^Ti{St+i) P{St+i\St) ■ ■■P{Si\So)^Tj{So) 

So,---,St+i 

= E i^T^iSt)i^TJ{So)P{St\St-l)■■■PiSl\So)MSoT (40) 

So,---,St 

Mc 

^ (2M,)-i E[V^T.(^.+*)^T,(5.) +^T^(^J^T,(5.+t)]. (41) 

CT=1 
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The first term in expression (^T]) follows immediately from expression (^OD; to obtain the 
second term one has to use the time reversal symmetry of a stochastic process that satisfies 
detailed balance, viz., 

PiSt+^\St) ■ ■ ■ P{S2\S,)MSi? = P{Si\S2) ■ ■ ■ PiSt\St^,)MSt+i)''. (42) 

Again, these estimators satisfy the zero-variance principle mentioned above, as long as 
the expressions are used as written, i.e., without symmetrization with respect to i and j. 



IV. TRIAL VECTORS 



As we mentioned, the form of the trial vectors used in these calculations is a major 
factor determining the statistical accuracy of the results. It not too difficult to make an 
initial guess for the form of the eigenvector corresponding to the second largest eigenvalue 
of the Markov matrix. Numerically exact calculations for small systems show that this 
eigenvector is antisymmetric under spin inversion, which is a manifestation of the longevity 
of fluctuations of the magnetization and not a peculiarity of small systems. 

This suggests the following initial approximation of the eigenvector belonging to the 
second largest eigenvector, the first sub-dominant eigenvector of the symmetrized Markov 
matrix P: 



ijTi{S) = mMS), (43) 

where m is the average magnetization. Figs. are a plots of ipi{S) / ^ipBiS) vs. the total 
magnetization M = L'^m for the exact eigenvector ipi computed for 3 x 3, 4 x 4, and 
5x5 nearest-neighbor Ising systems. For all three, the pre-factor m in Eq. (|4^) clearly 
captures a significant part of the truth, but there are two shortcomings. First of all, there is 
scatter, which indicates that ipi{S)/ipB{S) is a function of more than just the magnetization. 
Secondly, the "curve" is non-linear. The latter problem can be cured quite easily by replacing 
m in Eq. (^3[) by an odd polynomial m(l + a2m'^ + ■■■) with coefficients to be determined 
variationally. 

Similarly, computations for small systems (see Section |V A| for further details) suggest 



that the second largest eigenvalue is associated with an eigenvector that is even under spin 
inversion, as illustrated in Fig. ^. A trial vector of this form is readily constructed by 
replacing m on the right-hand side of Eq. (|43|) by a polynomial even in m. It turns out that 
the general picture as just described is largely independent of L. 

More in general, the plots shown in Figs. strongly suggest that the sub-dominant 
eigenvectors of the Markov matrix P, subject to the imposed spin, rotation and translation 
symmetries, are reasonably approximated by the Boltzmann distribution multiplied by a 
mode-dependent function of the magnetization. As can be seen in Figs. the number of 
nodes of this pre-factor increases by one as one steps down the spectrum, but it is also clear 
that, especially for the less dominant eigenvectors, the residual variance is significant. 

To begin to address the problem of the scatter and to improve the trial vector systemat- 
ically, it is necessary to identify other important variables besides the magnetization, and to 
incorporate them in the trial vector. We tried multi-spin correlations involving nearby spins 
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FIG. 1. Prc-factor V-'i/V'B of the first sub-dominant eigenvector of the Markov matrix vs. to- 
tal magnetization M for a 3 x 3 nearest-neighbor Ising {(3 = 0) lattice. For each value of M 
ipi{S)/ipB{S) is plotted for all configurations S with M{S) = M. All 512 points are plotted, but 
because of symmetries many coincide. 
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FIG. 2. Pre-factor t/'i/V'B of the first sub-dominant eigenvector of the Markov matrix vs. total 
magnetization M for a 4 x 4 nearest-neighbor Ising (/3 = 0) lattice. 
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FIG. 3. Pre-factor V-''i/'/''B of the first sub-dominant eigenvector of the Markov matrix vs. total 
magnetization M for a 5 x 5 nearest-neighbor Ising (/3 = 0) lattice. Although up to 400 configuration 
collapse onto a single point, crowding prevents individual resolution of most data points. 
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FIG. 4. Pre-factor i{j2/'ipB of the second sub-dominant eigenvector of the Markov matrix vs. total 
magnetization M for a 5 x 5 nearest-neighbor Ising (/3 = 0) lattice. 
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FIG. 5. Pre- factor ipa/i'B of the third sub-dominant eigenvector of the Markov matrix vs. total 
magnetization M for a 5 x 5 nearest-neighbor Ising (/3 = 0) lattice. 
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FIG. 6. Pre-factor ipi/ipB of the fourth sub-dominant eigenvector of the Markov matrix vs. total 
magnetization M for a 5 x 5 nearest-neighbor Ising (/? = 0) lattice. 
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FIG. 7. Pre-factor ip^/ipB of the fifth sub-dominant eigenvector of the Markov matrix vs. total 
magnetization M for a 5 x 5 nearest-neighbor Ising {/3 = 0) lattice. 

but after considerable failed experimentation we established that long-wavelength fluctua- 
tions of the magnetization are the suitable variables. This is reasonable when one compares 
e.g. the eigenvalue equations for ipo{S') and ipi{S') and realizes that the eigenvalues differ 
only very little from unity except for very small systems. We therefore used the Fourier 
components of the spin configuration, which are defined by 



where k = {ki, with < ki, ^2 < L; and Si-^i^ denotes the spin at lattice site (/i, I2). Note 
that m = mo,o- If we restrict ourselves to eigenvectors that are translationally invariant, the 
arguments presented in the previous paragraph yield the following trial odd or even vectors: 



The primes attached to the summation signs indicate that terms with kj = (0, 0) are ex- 
cluded. The coefficients a^, akik2' • • • polynomials in m, which are either odd or even 
under spin inversion and are to be chosen according to the desired symmetry. Rotation and 
reflection symmetries of the lattice are imposed by equating coefficients of the appropriate 
monomials in m. 

The results reported in Ref. || were obtained using a more complicated version of Eq. (^5D, 
namely 




(44) 



iPt{S) = iPb{S) [a^{m) + J2' a^^^^{m) my^^mk^ ^ki+ka.o 



ki,k2 




(45) 



ki,k2,k3 
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i)T{S) = ^b(5) [a (m) + flkikaM "^ki?^k2 4i+k2,o 

ki,k2 

+ J2' «kik2k3 ("^) "^kimkjmkg 5ki+k2+k3,0 + • • • 

ki,k2,k3 

X [a+(m) + Oj^^kaM "^ki^ka 5ki+k2,0 

ki,k2 

+ «kik2k3M"^kimk2mk3 4l+k2+k3,0 + • • •] (46) 

ki,k2,k3 

In those calculations also the coupling constant appearing in the Boltzmann factor was 
treated as variational parameter, but it turned out that the optimal value of this parameter 
was indistinguishable from the critical coupling. It does not seem that the more complicated 
form of expression (|46|) resulted in a major improvement, but we did not perform a systematic 
comparison of these trial vectors. 

The coefficients in the trial vector are treated as variational parameters. As in all non- 
linear fitting problems it is important to use parameters parsimoniously, and to do so one 
has to establish a hierarchy among these parameters. The scheme we used was to iterate 
the following step: (a) systematically add terms of increasing degree in m; (b) when this 
saturates increase the degree of terms with products of mk with mk 7^ (0, 0). 

The effectivity of this variational approach using low-momentum Fourier components, as 
described here, becomes apparent when one compares the variational eigenvalues with the 
exact numerical ones. For instance, the difference in the case of the second eigenvalue of the 
L = 5 nearest-neighbor model was only 2 x 10~^. 



V. NUMERICAL RESULTS 

A. Exact eigenvectors for small systems 

The full, symmetric Markov matrix P for an L x L Ising model is a 2^^ x 2^^ matrix, so 
that exact numerical calculations are possible only for very small systems; see e.g. results for 
L < 4 in Ref. 0. In the present work, we performed such exact computations for systems up 
to L = 5. In order to restrict the numerical task, we chose representations of P in subspaces 
with the appropriate symmetries. Two distinct symmetries were chosen, both of which 
impose invariance of the eigenvectors of P with respect to geometric translation, rotation 
and mirror inversion. The vectors were chosen to be either even or odd under spin inversion. 
This reduced by almost a factor 400 the dimensionality of P. In this way, the computation 
of a restricted set of eigenvectors became feasible for the resulting matrices of order 86,056 
for the L = 5 cases. For the diagonalization we made use of sparse-matrix methods and 



the conjugate- gradient method (see e.g., Ref. 19,20) which computes the eigenvector with 



the largest eigenvalue. Subsequent orthogonalization with respect to this eigenvector yields 
the eigenvector with the second largest eigenvalue, and further eigenvectors can be obtained 
similarly. Thus we obtained exact numerical solutions for six eigenvalues Xn and their 
corresponding eigenvectors ipi{S) (i = 0, ■ ■ ■ , 5) of the eigenvalue equation 

J2P{S,S')US') = Xu^l^,{S). (47) 

S' 
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The largest eigenvalue Alo is equal to 1, in accordance with the conservation of probability; 
its corresponding eigenvector satisfies i'oiS) = ipBiS), as follows from detailed balance. It 
is even with respect to spin inversion: ipo{S) = ipo{~S) where —S is obtained from S by 
inverting all spins. For all system sizes and models included here, we observed that the six 
leading eigenvectors, ordered according to magnitude of their eigenvalues, alternate between 
the odd and even subspaces: the first eigenvector is even, the second one is odd, the third 
one is even subspace, and so on, with the caveat that for L = 2 e.g. the odd subspace 
contains only two independent states. As we discussed above, the resulting eigenvectors 
provide useful information on how to construct trial vectors; moreover, the knowledge of 
accurate eigenvalues for L < 5 provided an powerful test of the Monte Carlo method, the 
results of which are presented in the following subsection. 



B. Monte Carlo calculations 



All simulations took place at the respective critical points of the models considered. 
This point is known exactly in the case of the nearest-neig hbor model [K^ = ln(l + V2)/2], 
and was determined numericallylll for the other two models: Kc = 0.1901926807 (2) for the 
equivalent-neighbor model, and = 0.6972207 (2) for the model with antiferromagnetic 



next-nearest neighbor interactions. The finite-size scaling analysis presented in Ref. p!2 
showed that, to the extent they are compatible with the numerical results, deviations from 
Ising universal behavior are extremely small. The raw simulation data used in this current 
paper include the data on which were based the numerical results for the largest relaxation 
time of the nearest-neighbor model, reported in Ref. ^. The latter results were obtained from 
8 X 10^ Monte Carlo samples for systems with finite sizes up to L = 15. The trial vector 
used for these computations consisted was of the form given in expression (^) and used 
up to 36 variational parameters. Also included in the present analysis are the simulations 
reported in Ref. |^, which contained 1.2 x 10^ Monte Carlo samples for all three models with 
system sizes up to L = 20. 

In addition, new simulations for each of the three models were performed, with a length 
of 2 X 10® Monte Carlo samples for system sizes up to L = 20, and of 1.6 x 10® Monte Carlo 
samples for system size L = 21. These new simulations used up to 89 variational parameters 
in the trial functions for each eigenvector of each model. 

In order to suppress biases due to deviations of randomness, we made use of a random 
number generator which combines two different binary shift registers such as described and 
discussed in Ref. 

The required Fourier components of the spatial magnetization distribution were sampled 
at intervals of one sweep for the smallest systems up to about 15 sweeps for the largest ones. 
The Monte Carlo calculation of the autocorrelation times (actually the eigenvalues of the 
Markov matrix) was performed for each run as a whole as well as separately for a number 
of up to 1024 blocks into which the run was split. This blocking procedure enabled us to 
estimate the statistical errors. Furthermore, the calculation of the eigenvalues according to 
V{t)N{t)~^ = A*^*) still depends on the time displacements t [see Eqs. (|39| ) and The 
calculation of A*^*^ was performed for time displacements tL^'^ = 0, 1,2, ■ ■ ■ up to 10 or 20 
of the above-mentioned intervals. For small t these eigenvalue estimates reflect variational 
bias due to the residual contributions of relaxation modes decaying faster than the mode 
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for which the trial vector was constructed. If the relaxation times of these faster modes are 
considerably shorter than that of the mode under investigation, one can clearly see a fast 
convergence of the eigenvalue estimate as a function of t. Convergence, however, occurs to 
a level that is only approximately constant because of the correlated statistical noise whose 
effect still depends on t. With increasing t, one can also observe that the statistical errors 
increase. The latter effect, which is as slow as the pertinent relaxation mode, occurs when the 
autocorrelations of the Monte Carlo sample are decreasing significantly with t. This situation 
was indeed observed for the largest eigenvalues; the data converged well with t before the 
coherence of the sampled data was lost. It was thus rather simple to select a "best estimate" 
of those eigenvalues. However, the situation for the smaller eigenvalues investigated here 
was much more difficult, because the relative differences between subsequent autocorrelation 
times are much smaller. The numerical results for the eigenvalues are listed in Appendix A. 



C. Determination of the dynamic exponent 

In two-dimensional Ising models, finite-size corrections are known that decay with finite 
size as L~^, and integral powers thereof may also be expected. In the absence of information 
on possible additional finite-size corrections of a different type that could occur in dynamic 
phenomena, we try to describe the finite-size data for the various autocorrelation times, as 



given in Eq. (pUf ), by the formula 



TL^L^Yl ^^L-'K (48) 

fc=0 

Here z is the dynamic exponent, the finite-size amplitude, and Uc is the number of 
correction-to-scaling terms included. Not explicitly shown in this notation is that the auto- 
correlation times depend on the relaxation mode and the model. 

On the basis of Eq. (^5]), a considerable number of least-squares fits were applied to 
the numerical results for the autocorrelation times. For each model and relaxation mode, 
one may vary both rzc, the number of correction terms, and the low-L cutoff specifying the 
minimal system size included in the fit. The smaller the number of corrections, the larger 
the low-L cutoff must be chosen in order to obtain an acceptable squared residual x^- A 
selection of fits that display the numerical trends is presented in Appendix B. 

The "best fits" were chosen on the basis of the criterion, the dependence on the low-L 
cutoff, and the mutual consistency of fits with different 71^.. The fits are summarized in Table 
|. Since the errors are not only of a statistical nature, but also depend on residual bias in the 
autocorrelation times and subjective choices made in the selection of the best fits, we quote 
error bars equal to two standard deviations as obtained from statistical considerations only. 
We believe that these two-sigma error estimates are conservative in the case of the analysis 
of the second and third largest eigenvalues of the /3 > models. The j3 = —1/4 model was 
found to be numerically less well-behaved: the statistical errors, as well as the corrections 
to scaling appear to be larger. Also, the construction of trial vectors was somewhat less 
successful than in the cases of the (3 > models. 

The new data are somewhat more accurate than, and consistent with our previous work.l3 
They are also consistent with the results of Wang and Hui for the slowest relaxation mode 
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TABLE I. Best estimates for the dynamic exponent z for five relaxation modes in three Ising-hke 
models. These results were selected from a much larger set of least-squares fits, obtained for 
different choices of the minimum system size and of the number of corrections taken into account 
(see Appendix B). The error estimate in the last decimal place of each entry is listed in parentheses, 
and is taken to be two standard deviations in the best fit. 



mode 


/3 = -l/4 


/3 = 


/3=1 


1 


2.164 


(3) 


2.1660 


(10) 


2.1667 


(5) 


2 


2.166 


(3) 


2.167 


(1) 


2.167 


(1) 


3 


2.164 


(5) 


2.170 


(2) 


2.167 


(1) 


4 


2.17 


(1) 


2.162 


(4) 


2.170 


(8) 


5 


2.15 


(2) 


2.17 


(1) 


2.17 


(1) 



of a different set of Ising-like models. They provide a clear confirmation of universality of the 
dynamic exponent, with regard to relaxation modes as well as models. Our best estimate 
z = 2.1667(5) applies to the slowest (odd) relaxation mode of the equivalent-neighbor 
{(3 = 1) model. 

This result for z is consistent with most of the recently published values. This agreement 
includes the results of Stauffeii on damage spreading in the Ising model. (The value listed 
in Ref. |^ was incorrectly quoted, for which we apologize.) It is slightly larger than the value 
2.14 derived by Alexandrowiczo on the basis of a scaling argument. 

Next we address the question whether the finite-size divergence of the autocorrelation 
times at criticality can be described by a dynamic exponent z = 2 when a logarithmic factor 
is included. This possibility was suggested by Domany,@ and pursued by SwendsenS and 
StaufferE3, who used very large lattices and found that this possibility seems inconsistent 
with one way of simulation, but not with a different one. Further references concerning this 
question are given in Ref. ^ 

Although the present work is restricted to very small system sizes, the data are relatively 
accurate. Thus we tried to fit the following form to the finite-size data for the slowest 
relaxation mode 

TL ^ L\l + MnL)(^ auL-^^). (49) 

fc=0 

Fixing z = 2 and taking 6 as a variable parameter, we found that this form could well 
describe the data for large enough n^. The quality, as determined by the criterion, of 
a number of such fits is shown in Table |T|. For comparison we include fits in which 2; is a 
variable parameter without a logarithmic correction. 

The results in Table ^ indicate that the fits with a variable exponent are usually better 
than those which include a logarithmic term, i.e., the residual decreases faster when the 
low-L cutoff is increased. This is especially apparent for small Uc and for the /3 = and 
(3 = 1 models, where the statistical accuracy is optimal. 

Finally we tried a fit according to Eq. (^) with both z and h as free parameters. The 
resolution of both parameters simultaneously is quite hard, and lies near the limit of what 
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TABLE II. Comparison between fits to the autocorrelation times, with and without a loga- 
rithmic finite-size dependence. These fits apply to the slowest relaxation mode. The first column 
shows the minimum system size included, the second the number of correction terms included. The 
fourth column displays the squared residual xl obtained when the dynamic exponent z was left 
free and the amplitude b of the logarithm was fixed at zero. The fifth column shows the squared 
residual xl when z was fixed at value 2, while b was left free. The sixth column lists the number 
of degrees of freedom of the fit for comparison. 



L > 


model 


Tic 


xl 


xl 


df 


5 


l3 = 


1 


239. 


274. 


76 


6 


/3 = 


1 


98. 


187. 


72 


7 


p = o 


1 


83.5 


127. 


67 


8 


(3 = 


1 


68.7 


87.9 


62 


9 


p = o 


1 


63.8 


71.9 


57 


10 


p = o 


1 


55.9 


57.5 


52 


8 


13 = 1 


1 


56.0 


71.1 


51 


9 


13 = 1 


1 


50.1 


57.4 


47 


10 


13 = 1 


1 


49.4 


49.7 


43 


4 


/? = 


2 


89.2 


341. 


76 


5 


13 = 


2 


85.0 


136. 


75 


6 


(3 = 


2 


83.3 


97.2 


71 



can be gleaned from the present data. For the nearest-neighbor model we find, using system 
sizes L = 8 and larger, and = 1, that z = 2.13 ± 0.07 and b = 0.05 ± 0.09, which again 
fails to support the presence of a logarithmic term. 



D. Universality of finite-size amplitudes 

In order to determine the leading amplitudes ao more precisely, we repeated the fits as 
used for the determination of the dynamic exponent, but with the value of the latter fixed 
at 2 = 13/6. We note that the combined results for z are consistent with this fraction. A 
considerable number of fits were made, and "best estimates" of the amplitudes are presented 



in Table HI. 



As mentioned in Ref. |^, according to a modest generalization of accepted ideas on uni- 
versality, the finite-size amplitudes of the autocorrelation times should satisfy ao = Aim^,, 
where the are non- universal, model-dependent constants; the subscript k refers to the 
specific model. We use the notation k = 2 + sgn(/3) so that k, = 1 refers to the model 
with ferromagnetic nearest- and antiferromagnetic next nearest neighbor couplings, k = 2 
denotes the nearest-neighbor model, and k = 3 to the equivalent-neighbor model. The Ai, 
i = 1, ■ ■ ■ ,5 are mode- dependent constants, whose ratios are universal. Since only the prod- 
uct matters, we are free to choose an arbitrary value for one of these constants. We chose 
to fix Ai = 1, so that all Ai should become universal constants. 
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TABLE III. Best estimates for the finite-size amplitudes of five relaxation modes in three 
Ising-like models. These results were selected from a much larger set of least-squares fits obtained 
for different choices of the minimum system size and of the number of corrections taken into account 
(see Appendix B). The error estimate in the last decimal place of each entry is listed in parentheses, 
and is taken to be two standard deviations of the best fit. The amplitudes can be written as the 
product of mode-dependent and model-dependent constants (see text); the difference in the last 
decimal place between the amplitudes and this product is shown between square brackets. 



mode/model 


K = l 


K = 2 


K = 3 


1 (odd) 


6.763 


(6) 


[-10] 


4.4089 


(13) 


[8] 


2.8312 


(5) 


[0] 


2 (even) 


0.2516 


(2) 


[2] 


0.16364 


(5) 


[0] 


0.10510 


(2) 


[0] 


3 (odd) 


0.1188 


(1) 


[1] 


0.07727 


(3) 


[0] 


0.04963 


(1) 


[0] 


4 (even) 


0.07195 


(7) 


[3] 


0.04677 


(4) 


[-4] 


0.03008 


(3) 


[2] 


5 (odd) 


0.0466 


(1) 


[-1] 


0.03041 


(3) 


[-1] 


0.01956 


(2) 


[2] 



The remaining constants (i = 2, ■ • ■ , 5) and were fitted accordingly to the ampli- 
tudes listed in Table |T|. The result of this least-squares fit is mi = 6.773 ± 0.003, m2 = 
4.4081 ±0.0009, mg = 2.8312±0.0005, A2 = 0.037123 ±0.000008, A3 = 0.017530±0.000004, 

= 0.010618 ± 0.000006, and A^ = 0.006901 ± 0.000005. This fit has 8 degrees of freedom 
and = S-Oj ^ good agreement with the assumptions of dynamic universality, and sug- 
gesting that our two-sigma error estimates are not unrealistic. The differences between our 
amplitude estimates and the fitted product AifriK are included in Table fT|, in units of the 
last decimal place listed. 



VI. DISCUSSION 

The analysis presented above produces an apparently highly accurate estimate of the 
dynamic critical exponent z = 2.1667±0.0005 for the case of the equivalent neighbor model, 
where the statistical accuracy and the convergence are best. However, even in this case it 
is not possible to rule out a divergence of the relaxation time of the form L^(l ± blnL). 
Nevertheless, our results viewed in their totality make this behavior rather unlikely. Firstly, 
the assumption of this logarithmic form yields fits that converge less rapidly, as mentioned 
above. Secondly, one would have to have 6 1/6 universally, independent of model and 
relaxation mode, since our results for the range of system size we studied are consistent with 
a divergence of the form ti (x (x L^^/^ ^ L^(l ± 1/6 In L). Universality of the amplitude 
of a logarithmic correction would be quite unusual and does not fit into any theoretical 
framework of which we are aware. 

Some open questions remain regarding the optimized variational vectors to which this 
computation owes its accuracy. Variational basis vectors of the general form given in Eq. (^5]) 
are special in the sense that all parameters enter linearly. The method outlined here does 
not require this feature and in fact it was not present in previous computations, reported in 
Ref. |. 
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For most of the results presented here we used trial vectors with linear parameters opti- 
mized by minimization of the variance of the configurational eigenvalues. As an apparently 
equivalent alternative, the full set of symmetrized monomials in the Fourier coefficients of 
the spin configuration could be chosen as basis vectors rather than the linear combinations 
defined in Eq. (|i5|). With this choice, the basis vectors would not have contained any pa- 
rameters, but employing this bigger truncated basis, the same linear parameters would have 
been reintroduced by computing the matrix elements Nij and Vij and solving the generalized 



eigenvalue problem defined by Eq. (pSj). In this way, we would have obtained the coefficients 
for which the Rayleigh quotient is stationary, at least if the summation over configurations 
could have been done exactly. Proceeding in this way, we could have altogether skipped the 
optimization scheme based on minimization of the variance of the configurational eigenvalues 
[cf. Eq. m) and Section |irB|] . 



The obvious question is what is accomplished by the non-linear minimization of the 
variance of the configurational eigenvalues. We do not yet have a convincing answer to this 
question. On the one hand, it is not difficult to show that the zero variance principle holds 
for individual eigenstates. That is, if an eigenvector can represented exactly as a linear 
combination of the basis vectors, the variational (t = 0) case will already produce the exact 
result even if it the eigenvalue problem contains eigenvectors that cannot be represented 
exactly in the truncated basis. This is in fact precisely what happens for the dominant even 
eigenvector: this vector is represented exactly even in the truncated basis we use, and indeed 
its eigenvalue is reproduced exactly. On the other hand, our tentative numerical experiments 
show that the optimization method produces more accurate results, which is not surprising 
when one considers that the optimized basis functions give rise to much smaller truncated 
basis sets. This in turn yields a generalized eigenvalue problem involving much smaller 
matrices that are numerically and statistically much more robust. 

A related problem with a large basis set is that the matrix N in the generalized eigenvalue 
problem of Eq. (^) becomes numerically singular. In fact, the only way in which we were 
able to obtain meaningful results at all is by performing the inversion in the usual regularized 
fashion as follows. Use the fact that N is symmetric and non-negative definite to write it in 
the form 

N = Wdmg{fil,...,fil)Wl (50) 
Then define a regularized inverse of iV^ as follows 

iv-^ = ^^diag(/i^^...,/x;^)w^^ (51) 

where fi~^ = if /ij exceeds a suitable chosen threshold, e.g. the square root of the machine 
accuracy, and fi'^^ = otherwise. The non- vanishing eigenvalues of N~2PN^2 then yield a 
subset of the eigenvalues of P that are least affected by the numerical singularity of N. 

Although further modifications of the computational procedures may lead to additional 
improvements of our technique, the numerical results obtained thus far are already quite 
promising and the question arises what further applications are obvious in the field of dy- 
namics of Monte Carlo methods. 

For instance, it seems well possible to apoly the present techniques in three dimensions, 
and to spin-conserving Kawasaki dynamicsEI although it is clear that the construction of 
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trial vectors will have to be modified. In the same context, we note that direct application 
of the method used in this paper to the dynamics of cluster algorithms^ is frustrated by 
the requirement that one should be able to compute {S\P\iPt) numerically exactly. An 
additional problem for such dynamics from the perspective of our approach is that the 
concept correlation time has to be handled carefully in this context. 

Let us demonstrate this point by means of the following thought experiment: the appli- 
cation of the Wolff algorithm to the ferromagnetic, critical Ising model. As usual, we define 
the autocorrelation times in terms of the eigenvalues of the stochastic matrix. In order to 
enable a comparison with other types of dynamics, we choose our unit of time as L'^'^'^y^ 
Wolff steps {L is the linear system size, d the dimensionality and yh the magnetic renormal- 
ization exponent). Since the average Wolff cluster consists of a number of sites proportional 
to L'^yh-'^^ this choice guarantees that, under equilibrium conditions, an average number of 
order L'' spins is processed per unit of time. 

Because of the efficiency of the Wolff algorithm, only a few units of time are needed to 
generate an independent spin configuration under practical circumstances. However, if the 
fully ordered antiferromagnetic state is chosen as the initial spin configuration, a number of 
Wolff steps of order L*^ is required to remove the the antiferromagnetic order, i.e., its relax- 
ation to equilibrium is anomalously slow. A less extreme but related phenomenon is observed 
under practical Wolff simulation conditions at equilibrium: from time to time large critical 
fluctuations occur that bring the system in a state of relatively large disorder and small 
magnetization. These configurations are relatively long-lived. In the time-autocorrelation 
function of the magnetization, this phenomenon translates into a slower-than-exponential 
decayed, at least on the numerically accessible time scale. In the language of Eq. ( |T0| ) such 
a situation follows if one assumes the existence of anomalously large autocorrelation times 
Tii associated with anomalously small amplitudes q. Under these circumstances we cannot 
exclude the possibility that the longest relaxation time following from the Markov matrix 
for the Wolff simulation of a finite system corresponds with an extremely unlikely deviation 
from equilibrium. Since this kind of fluctuations may have too low a probability to be of 
practical significance, these considerations suggest the possibility that the time needed to 
generate an 'independent configuration' is not simply related to the second largest eigenvalue 
of the Markov matrix, but rather to some intricate average, possibly involving the complete 
spectrum. 
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APPENDIX A: NUMERICAL RESULTS FOR THE AUTOCORRELATION 

TIMES 



This Appendix contains Table [TVl with numerical results for the eigenvalues of the Markov 
matrix for five relaxation modes of three Ising-like models as described in Section ^ In order 
to reduce the amount of data, results of different runs for the same model and system size 
are combined into a single entry. The autocorrelation times are directly related to these 
data via Eq. 
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APPENDIX B: NUMERICAL FITS INVOLVING Z 



This Appendix lists Table |V] containing a selection of fits used to determine the dynamic 
exponent z. These fits were made for five relaxation modes of three Ising-like models. For 
a given number of corrections-to-scaling n^. the minimum system size must not be too small 
in order to get an acceptable fit as becomes apparent on the basis of the criterion. The 
Table contains a number of entries illustrating this point. It includes data for the amplitudes 
ctQ which support the hypothesis of universal amplitude ratios. However, statistically more 
accurate amplitudes (see Table HD can be determined on the basis of the assumption z = 
13/6 which is consistent with the results in Table |V|. 
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TABLE IV. Eigenvalues of the stochastic matrix for different relaxation modes and Ising-Iike 
models. The first column indicates the relaxation mode, the second one refers to the linear 
system size. Next follow the eigenvalues for the three models, identified by the sign of the 
next-nearest-neighbor coupling determined by k. Each eigenvalue is followed by the statistical 
error in its last decimal place between parentheses. Also shown are exact numerical results for 
system sizes up to L = 5. The latter entries are thought to be accurate in all decimal places given. 
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TABLE V. Fits used for the determination of the dynamic exponent z. These fits were made 
for different relaxation modes (i, column 1) of three distinct models identified by their value of k 
(column 2). The third column indicates the minimum system size, and the fourth one the number of 
correction terms rtc that were included in the fit. The following columns list the estimated dynamic 
exponent, its statistical uncertainty (two standard deviations), the residual x^, the number of 
degrees of freedom, and the amplitude of the term describing the pertinent relaxation mode. 
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